%% Contains all the plotting functions for easy acces
% Plotmode is used to select which plot to perform

function varargout = GeneralPlotting(plotmode,savevars,cfg,exvars,depvars,indepvars,r,t)
%% Creates the plot for each region
global globals % Used to keep track of earlier iteration and for realtime plots

modes=indepvars.UseModes(r,:);
if nargout > 0
        if ~isfield(cfg,'Run')
            varargout{1}     = [ 'test ' datestr(now,'yyyy-mm-dd HH-MM-SS' ) ' R ' cfg.allregions{r}];
        else
            varargout{1}     = [ 'test ' datestr(now,'yyyy-mm-dd HH-MM-SS' ) ' R ' cfg.allregions{r} ' Run ' num2str(cfg.Run)];
        end
end

switch plotmode
    case {'PF or SW over time'}       
        if ishandle(4);set(4,'Units','centimeter'); cornerpos = get(4,'Position'); close(4);end;
        cfg.savename        = ['test ' datestr(now,'yyyy-mm-dd HH-MM-SS' ) ' R ' cfg.allregions{r}];
        [bs1 bs2]           = mkdir(cfg.outputfolder);

        globals.fid = figure(4);set(4,'Units','centimeter')
        PlotRange   = 1:130;
        if ~exist('cornerpos','var'); cornerpos = get(globals.fid,'Position'); end
        
        if strcmp(cfg.Model,'ModelGcam')
            cfg.WeightFields = 'ShareWeight';
            cfg.Power        = -2;
            plot(squeeze((depvars.(cfg.WeightFields)(r,:,3))./squeeze(depvars.(cfg.WeightFields)(r,:,4))).^cfg.Power)

        else
            cfg.WeightFields = 'PreferenceFactor';
            cfg.Power        = 1;
            plot(squeeze(depvars.(cfg.WeightFields)(r,:,:)).^cfg.Power)
        end
        title(['Development toward saturation ' cfg.allregions{r}]);
        ylabel([cfg.WeightFields ' to the power of ' num2str(cfg.Power)]);
        set(gca,'XLim',[1 130],'XTick',PlotRange(round(linspace(PlotRange(1),PlotRange(end),4))),'XTickLabel',num2cell(cfg.years(PlotRange(round(linspace(PlotRange(1),PlotRange(end),4))))))
    
        
        
    case {'Barplot errors'}
       if ishandle(4);set(4,'Units','centimeter'); cornerpos = get(4,'Position'); close(4);end;
        cfg.savename        = ['test ' datestr(now,'yyyy-mm-dd HH-MM-SS' ) ' R ' cfg.allregions{r}];
        [bs1 bs2]           = mkdir(cfg.outputfolder);

        if r == cfg.CalculatedRegions(1); globals.barcontent = zeros([12 1]); end

        globals.fid        = figure(4);set(4,'Units','centimeter')
        selectedfields{1}     = 'ModeFit'; n=1;

        globals.barcontent(r) = depvars.errors.(selectedfields{n})(r);
        cfg.shortregions = {'NAM','PAO','WEU','EEU','FSU','Afr','CPA','LAM','MEA','PAS','SA','Wor'};
        globals.fid = barh(globals.barcontent);

        set(gca,'yticklabel',cfg.shortregions)
        xlabel([selectedfields{n}]); ylabel('Region'); title(['Comparison for ' selectedfields{n}]);

        configs.savename = ['Barh ' selectedfields{n} ' ' datestr(now,'yyyy-mm-dd HH-MM-SS' )];
                
    case 'Single Plot'
        ErrorRange = 1:35;
        PlotRange  = 1:130;
        if ishandle(4);set(4,'Units','centimeter'); cornerpos = get(4,'Position'); close(4);end;
        cfg.savename        = ['test ' datestr(now,'yyyy-mm-dd HH-MM-SS' ) ' R ' cfg.allregions{r}];
        [bs1 bs2]           = mkdir(cfg.outputfolder);

        globals.fid = figure(4);set(4,'Units','centimeter')
        if ~exist('cornerpos','var'); cornerpos = get(globals.fid,'Position'); end
        plot(permute(exvars.TravelDemand(r,ErrorRange,:),[2 3 1]),'--');
        hold on    
        plot(permute(depvars.TravelDemand(r,PlotRange,:),[2 3 1]));
        hold off
        title(['TravelDemand ' cfg.allregions{r}]);
        ylabel('Demand in PKM/yr');
        legend(cfg.modes);
        set(gca,'XLim',[1 130],'XTick',PlotRange(round(linspace(PlotRange(1),PlotRange(end),4))),'XTickLabel',num2cell(cfg.years(PlotRange(round(linspace(PlotRange(1),PlotRange(end),4))))))
    
    case 'ProjectionArea'
        %% Replicates the plot from the GCAM publication
        
        if ishandle(4);set(4,'Units','centimeter'); cornerpos = get(4,'Position'); close(4);end;
        globals.fid = figure(4);set(4,'Units','centimeter')
        if ~exist('cornerpos','var'); cornerpos = get(globals.fid,'Position'); end
        PlotRange = 36:130; 

        area(squeeze(depvars.TravelDemand(r,PlotRange,[3 2 4])) .* repmat(exvars.Population(r,PlotRange)',[1 3])/10^6);
        child = get(gca,'Children'); childcolor= ([0 0 0;1 1 1;  0.6 0.6 0.6; 0.9 0.9 0.9;]);
        for n=1:3; set(child(end-n+1),'FaceColor',childcolor(n,:)); end
        pubgraph(globals.fid,10,1.5,'white') % makes the graphs a little more readable
        title(['TravelDemand ' cfg.allregions{r} ', Run:' num2str(cfg.Run) ', CD:' num2str(indepvars.CostDistrParameter(r),3) ' EY:' num2str(indepvars.PriceElasticity(r),3) ' IY:' num2str(indepvars.IncomeElasticity(r),3)]);
        title(cfg.allregions{r});
        set(gca,'XTick',cfg.Time(round(linspace(1,length(PlotRange),7))),'XTickLabel',num2cell(cfg.years(PlotRange(round(linspace(1,length(PlotRange),7))))))
        ylabel('TravelDemand total, trillion pkm');
        plotranges = [40 14 120 160];
        if r<=4
        set(gca,'YLim',[0 plotranges(r)]);
        end
%         h=legend({'LDV','PT','HSM'},'Location','NorthWest'); 

    case {'Regression and Scatter'}
        %% Replicates the plot from the GCAM publication
        [cfg,exvars,depvars,indepvars] = TransformModesandRegions(cfg,exvars,depvars,indepvars,'regions11');
        [cfg,exvars,depvars,indepvars] = TransformModesandRegions(cfg,exvars,depvars,indepvars,'modes4');    
        
        if ishandle(4);set(4,'Units','centimeter'); cornerpos = get(4,'Position'); close(4);end;
        globals.fid = figure(4);set(4,'Units','centimeter')
        if ~exist('cornerpos','var'); cornerpos = get(globals.fid,'Position'); end
        PlotRange = 36:130;
        varargout    = {['test ' datestr(now,'yyyy-mm-dd HH-MM-SS' ) ' R ' cfg.allregions{r}]}; 
        
        plotregression(reshape(exvars.TravelDemand(1,1:35,:),[1,4*size(exvars.TravelDemand,2)]),...
            reshape(depvars.TravelDemand(1,1:35,:),[1,4*size(exvars.TravelDemand,2)]));
        title(['TravelDemand ' cfg.allregions{r}]);
        set(gca,'XTick',cfg.Time(round(linspace(1,length(PlotRange),7))),'XTickLabel',num2cell(cfg.years(PlotRange(round(linspace(1,length(PlotRange),7))))))
        ylabel('TravelDemand total pkm');
        h=legend({'Bus+Train','Car','HS Train+Air'}); 
        pubgraph(globals.fid,10,1.5,'white') % makes the graphs a little more readable
        
    case {'TTB and TMB'}
        if ishandle(4);set(4,'Units','centimeter'); cornerpos = get(4,'Position'); close(4);end;
        globals.fid = figure(4);set(4,'Units','centimeter')
        if ~exist('cornerpos','var'); cornerpos = get(globals.fid,'Position'); end
        if ~isfield(cfg,'Run')
        cfg.savename        = [ 'test ' datestr(now,'yyyy-mm-dd HH-MM-SS' ) ' R ' cfg.allregions{r}];
        else
            cfg.savename    = [ 'test ' datestr(now,'yyyy-mm-dd HH-MM-SS' ) ' R ' cfg.allregions{r} 'Run:' num2str(cfg.Run)];
        end
        
        [bs1 bs2]           = mkdir(cfg.outputfolder);

        globals.fid = figure(4);set(4,'Units','centimeter')
        PlotRange = 1:130;
        ErrorRange= 1:35;
            
%         subplot(2,1,1);
        plot(permute(indepvars.TTB(r,PlotRange,:),[2 3 1]),'--','Color',[0.2 0.2 0.2]);
        hold on 
        plot(permute(sum(depvars.TravelDemand(r,PlotRange,:)./exvars.Speed(r,PlotRange,:),3),[2 3 1]),'Color',[0.2 0.2 0.2]);
%         plot(permute(sum(exvars.TravelDemand(r,1:35,:)./exvars.Speed(r,1:35,:),3),[2 3 1]),'r');
        hold off
        title(['TTB pred. vs model ' cfg.allregions{r} ' [Hours/year]']);
        set(gca,'XTick',PlotRange(round(linspace(PlotRange(1),PlotRange(end),7))),'XTickLabel',num2cell(cfg.years(PlotRange(round(linspace(PlotRange(1),PlotRange(end),7))))))
        ylabel('TTB');
%             
%         subplot(2,1,2);
%         plot(permute(sum(depvars.TravelDemand(r,PlotRange,:).*exvars.TravelPrice(r,PlotRange,:),3)./exvars.Income(r,PlotRange),[2 3 1]),'--','Color',[0.2 0.2 0.2]);
%         hold on 
%         plot(permute((exvars.TMBCorrection(r,1) + 0.035 + 0.085 .* sum(depvars.ModeSplit(r,PlotRange,cfg.ModesFasterThanLDV),3)),[2 3 1]),'Color',[0.2 0.2 0.2]);
%         hold off
%         title(['TMB pred. vs model ' cfg.allregions{r} ' [Hours/year]']);
%         set(gca,'XTick',PlotRange(round(linspace(PlotRange(1),PlotRange(end),7))),'XTickLabel',num2cell(cfg.years(PlotRange(round(linspace(PlotRange(1),PlotRange(end),7))))))
%         ylabel('TMB');

        pubgraph(globals.fid,10,1.5,'white') % makes the graphs a little more readable
        set(globals.fid,'Position',get(globals.fid,'PaperPosition').*[0 0 1 1]+cornerpos.*[1 1 0 0]);

        if isequal(savevars,'save')
        savestuff(cfg)
        end
        
    
    case {'HistoricPlot','ProjectionPlot','AllYearsPlot'}
        %% Line plots with all fitting information
            switch plotmode
                case 'HistoricPlot'
                PlotRange = 1:35;
                ErrorRange= 1:35;
                
                case 'ProjectionPlot'
                PlotRange = 36:130;
                ErrorRange= 1:35;

                case 'AllYearsPlot'
                PlotRange = 1:130;
                ErrorRange= 1:35;
            end
            
            if ishandle(4);set(4,'Units','centimeter'); cornerpos = get(4,'Position'); close(4);end;
            cfg.savename        = ['test ' datestr(now,'yyyy-mm-dd HH-MM-SS' ) ' R ' cfg.allregions{r}];
            [bs1 bs2]           = mkdir(cfg.outputfolder);

            globals.fid = figure(4);set(4,'Units','centimeter')
            if ~exist('cornerpos','var'); cornerpos = get(globals.fid,'Position'); end
subplot(3,2,1);
            plot(permute(exvars.TravelDemand(r,ErrorRange,:),[2 3 1]),'--');
            hold on    
            plot(permute(depvars.TravelDemand(r,PlotRange,:),[2 3 1]));
            hold off
            title(['TravelDemand ' cfg.allregions{r}]);
            ylabel('Demand in PKM/yr');
            set(gca,'XTick',PlotRange(round(linspace(PlotRange(1),PlotRange(end),7))),'XTickLabel',num2cell(cfg.years(PlotRange(round(linspace(PlotRange(1),PlotRange(end),7))))))
%             set(gca,'YScale','log')

subplot(3,2,2);
            tempTD = exvars.TravelDemand(r,ErrorRange,:);
            tempTD(1,:,~modes) = 0;
            exvars.ModeSplit = tempTD(:,:,:) ./ repmat(sum(tempTD(:,:,:),3),[1 1 length(cfg.modes)]);
            plot(permute(depvars.ModeSplit(r,PlotRange,:),[2 3 1]));
            hold on    
            plot(permute(exvars.ModeSplit,[2 3 1]),'--');
            title(['ModeSplit ' cfg.allregions{r}]);
            hold off
            ylabel('Fraction per mode');
            h=legend(cfg.modes); 
            set(gca,'XTick',PlotRange(round(linspace(PlotRange(1),PlotRange(end),7))),'XTickLabel',num2cell(cfg.years(PlotRange(round(linspace(PlotRange(1),PlotRange(end),7))))))
            set(h,'Position',[0.85 0.06 0.13 0.29])
            
subplot(3,2,3);
            switch cfg.Model
                case {'ModelTravelLinearTW','ModelTravelOptimizedTW'}
                plot(permute(depvars.RelativeModeCostTravel(r,PlotRange,:),[2 3 1]));
                title(['Relative ModeCost ' cfg.allregions{r} ' [% deviation from mean]']);
                ylabel('Price per mode');
                set(gca,'XTick',PlotRange(round(linspace(PlotRange(1),PlotRange(end),7))),'XTickLabel',num2cell(cfg.years(PlotRange(round(linspace(PlotRange(1),PlotRange(end),7))))))
            

                case 'ModelGcam'
                plot([permute(depvars.ModeCostGcam(r,PlotRange,:),[2 3 1])]);
                hold on
                plot(sum(depvars.ModeCostGcam(r,PlotRange,:).*depvars.ModeSplit(r,PlotRange,:),3),'--')
                hold off
                title(['ModeCost ' cfg.allregions{r} ' [USD/PKM]']);

                case { 'ModelPoles'}
                hold on 
                plot(depvars.EquipmentRate(r,PlotRange,3));
                plot(exvars.LoadFactor(r,PlotRange,3),'--');
                hold off
                title(['Car Equipment rate ' cfg.allregions{r} ' [v/p]']);
                ylabel('cars/p');
                set(gca,'XTick',PlotRange(round(linspace(PlotRange(1),PlotRange(end),7))),'XTickLabel',num2cell(cfg.years(PlotRange(round(linspace(PlotRange(1),PlotRange(end),7))))))
            

            end
            
            ylabel('Price per mode');
            set(gca,'XTick',PlotRange(round(linspace(PlotRange(1),PlotRange(end),7))),'XTickLabel',num2cell(cfg.years(PlotRange(round(linspace(PlotRange(1),PlotRange(end),7))))))
            
            if isfield(depvars,'TimeWeight') && depvars.TimeWeight (r,1) ~= depvars.TimeWeight(r,2);
            end
            
subplot(3,2,4);
          switch cfg.Model
                case {'ModelGcam'}
                plot(sum(exvars.TravelDemand(r,ErrorRange,indepvars.UseModes(r,:,:)),3),'--');

                hold on 
                plot(sum(depvars.TravelDemand(r,PlotRange,indepvars.UseModes(r,:,:)),3));

                hold off
                title(['TravelDemand actual vs normal ' cfg.allregions{r} ' [Hours/year], TD data vs world (right axis)']);
                ylabel('Price per mode');
                set(gca,'XTick',PlotRange(round(linspace(PlotRange(1),PlotRange(end),7))),'XTickLabel',num2cell(cfg.years(PlotRange(round(linspace(PlotRange(1),PlotRange(end),7))))))
            
                case {'ModelPoles','ModelSimple'}
                    plot(permute(indepvars.TTB(r,PlotRange,:),[2 3 1]),'--','Color',[0.2 0.2 0.2]);
                    hold on 
                    plot(permute(sum(depvars.TravelDemand(r,PlotRange,:)./exvars.Speed(r,PlotRange,:),3),[2 3 1]),'Color',[0.2 0.2 0.2]);
                    hold off
                    title(['TTB pred. vs model ' cfg.allregions{r} ' [Hours/year]']);
                    set(gca,'XTick',PlotRange(round(linspace(PlotRange(1),PlotRange(end),7))),'XTickLabel',num2cell(cfg.years(PlotRange(round(linspace(PlotRange(1),PlotRange(end),7))))))
                    ylabel('TTB');
                
                case {'ModelTravelLinearTW','ModelTravelOptimizedTW'}
                    plot(depvars.TimeWeight(r,PlotRange));
                    title(['TimeWeight ' cfg.allregions{r}]);
                    ylabel('TimeWeight');
                    set(gca,'XTick',PlotRange(round(linspace(PlotRange(1),PlotRange(end),7))),'XTickLabel',num2cell(cfg.years(PlotRange(round(linspace(PlotRange(1),PlotRange(end),7))))))
          end
            if isfield(cfg,'CalibrationArgs')
subplot(3,2,5);
            cfg.axiscontent = {};
            cfg.barcontent  = [];
            calargs         = unique(cfg.CalibrationArgs{r});
            for n = 1:length(calargs);
                eval(['elements = size(' calargs{n} ',3);'])
                if elements == length(cfg.modes)
                    for m= 1:length(cfg.modes)
                        if eval([calargs{n} ' (r,1,m) ~=0' ]) 
                        cfg.axiscontent{end+1} = strrep(strrep(strrep(strrep([calargs{n} ' ' cfg.modes{m}],'indepvars.',''),'depvars.',''),'PreferenceFactor','PF'),'ShareWeight','SW');
                        eval(['cfg.barcontent(1,end+1)  = ' calargs{n} '(r,1,m);']);
                        end
                    end
                else
                    cfg.axiscontent{end+1} = strrep(strrep(strrep(strrep([calargs{n}],'indepvars.',''),'depvars.',''),'PreferenceFactor','PF'),'ShareWeight','SW');
                    eval(['cfg.barcontent(1,end+1)  = ' calargs{n} '(r,1,1);']);
                end
            end
            barh(1:length(cfg.barcontent),cfg.barcontent,'r');
            set(gca,'YTickLabel',cfg.axiscontent); 
            title('Calibration parameters');
            ylim([0.5 length(cfg.barcontent)+0.5]);
            end
            
            cfg.h.textbox = annotation('textbox',[ 0.52 0.02 0.30 0.29 ], 'String','hi','HorizontalAlignment','center','FontSize',8);
            set(cfg.h.textbox,'String', [datestr(now,'yyyy-mm-dd HH-MM-SS'),['Model: ' cfg.Model], ['Solver: ' cfg.Solver], ['NRSME OF/TTB/MODE/TD: '] , [ num2str(depvars.errors.OFfit(r),2) '/' num2str(depvars.errors.TTBfit(r),2) '/' num2str(depvars.errors.ModeFit(r),2) '/' num2str(depvars.errors.TDfit(r),2)]; cellstr(cfg.plotinfo{r})],'EdgeColor',[1 1 1])
            pubgraph(globals.fid,10,1.5,'white') % makes the graphs a little more readable
            set(globals.fid,'Position',get(globals.fid,'PaperPosition').*[0 0 1 1]+cornerpos.*[1 1 0 0]);            
            
            if isequal(savevars,'save')
            savestuff(cfg)
            end
               
        
    case 'SaveToGlobal'
        %% Allows cetrain veriables to be saved to the global workspace
        if isfield(cfg,'CalibrationArgs')
        globals.Error(end+1)                  = mean(depvars.Resid(r,t));
        for n=1:length(cfg.CalibrationArgs{r,1})
        eval(['globals.CalibrationFactors.' (cfg.CalibrationArgs{r,1}{n}) '(:,end+1)' ...
        '=' cfg.CalibrationArgs{r,1}{n} '(r,1,:);']);
        end
        end
   
    case 'PlotMinimize'
        %% PlotMinimize
        if cfg.PlotMinimize==1; globals.fid=figure(1); 
            subplot(3,2,1); plot(real(globals.Error)); title([cfg.allregions{r} ', Time = ' num2str(t(end)) ', Residual^2 = ' num2str(depvars.Resid(r,t(end)))]);
            subplot(3,2,2); bar(squeeze(depvars.ModeSplit(r,t(end),:))); title('ModeSplit')
            
            subplot(3,2,3); bar([squeeze(exvars.TravelDemand(r,t(end),:)) squeeze(depvars.TravelDemand(r,t(end),:))]); title('TravelDemand');
            for n=1:length(cfg.CalibrationArgs{r,1})
            subplot(3,2,3+n); 
            eval(['plot(permute(globals.CalibrationFactors.' cfg.CalibrationArgs{r,1}{n} ',[2,1]));']); title(cfg.CalibrationArgs{r,1}{n}); legend([cfg.modes]); legend('boxoff');
            end

            pubgraph(globals.fid,10,2,'white') % makes the graphs a little more readable
            set(globals.fid,'Units','centimeter')
            set(globals.fid,'Position',get(globals.fid,'PaperPosition').*[0 0 1 1]+get(globals.fid,'Position').*[1 1 0 0]);
            disp('Any key to continue'); pause;
        end
    case 'LiveMinimizeVisuals'
        
        GeneralPlotting('SaveToGlobal','',cfg,exvars,depvars,indepvars,r,t)
        plotcols = ceil((length(cfg.CalibrationArgs{r,1})+3)/2);
        
        globals.fid = figure(2);
        set(globals.fid,'Color','w');
        subplot(plotcols,2,1);
        plot(squeeze(real(depvars.TravelDemand(r,t,:))));
        title('TravelDemand exvar vs depvar')
        set(gca,'YScale','log');

        subplot(plotcols,2,2);
        plot(real(squeeze((depvars.TravelDemand(r,t,:) - exvars.TravelDemand(r,t,:)).^2)));
        title('Squared difference TravelDemand Schafer and model'); 


        subplot(plotcols,2,3);
        plot(real(globals.Error));
        title('Total error development');

        for n=1:length(cfg.CalibrationArgs{r,1})
            subplot(plotcols,2,3+n);
            eval(['plot(permute(globals.CalibrationFactors.' cfg.CalibrationArgs{r,1}{n} ',[2,1]));'])
            title(cfg.CalibrationArgs{r,1}{n})  
        end
        set(globals.fid,'Units','centimeter')
        set(globals.fid,'Position',get(globals.fid,'PaperPosition').*[0 0 1 1]+get(globals.fid,'Position').*[1 1 0 0]);
        
    case 'GUIPlots'
        %% GUIPlots for the Graphical user interface
        if cfg.UpdateGUIPlot == 1
            figure(3); subplot(2,2,2); plot([depvars.Resid(r,1:1:(t<36)*(t-35)+35,:) zeros([1,(t>35)*(t-35)])]);title('Residual error')
            subplot(2,2,4); plot([permute(exvars.TravelDemand(r,1:(t<36)*(t-35)+35,:),[2 3 1]); zeros([(t>35)*(t-35),7])],'--'); 
            hold on; 
            plot(permute(depvars.TravelDemand(r,1:t,:),[2 3 1])); title('TravelDemand per mode')
            hold off; 
            set(gca,'YScale','log');
            pause(0.05);
        end
        
    case 'Send to Workspace'

        assignin('base', 'cfg', cfg);
        assignin('base', 'exvars', exvars);
        assignin('base', 'depvars', depvars);
        assignin('base', 'indepvars', indepvars);
            
    case 'Sensitivity Analysis'
        %% Makes a 3d-plot of two parameters after optimization
        
        xfields                 = fieldnames(indepvars);
        selection               = SelectionMenu(fieldnames(indepvars),'Select error indicators to compare');
        cfg.XFields             = xfields(selection);
        
        z=1;
        for n=1:length(cfg.XFields)
            if size(indepvars.(cfg.XFields{n}),3)~=1 && sum(indepvars.(cfg.XFields{n})(r,1,:)~=indepvars.(cfg.XFields{n})(r,1,1))~=1;
                loopmodes = find(indepvars.UseModes(r,:));
                usemodes  = cfg.modes;
                spacer    = ' ';
            elseif sum(indepvars.(cfg.XFields{n})(r,1,:)~=indepvars.(cfg.XFields{n})(r,1,1))==1
                loopmodes = find(indepvars.(cfg.XFields{n})(r,1,:)~=indepvars.(cfg.XFields{n})(r,1,1));
                usemodes  = cfg.modes;
                spacer    = ' ';
            else
                loopmodes = 1;
                usemodes  = {''};
                spacer    = '';
            end
            for m=loopmodes
                varyvars                        = indepvars;
                datapoints                      = [40];
                percentualrange                 = linspace(0.025,2,datapoints(1));
                output                          = {};
                varyvars.(cfg.XFields{n})       = indepvars.(cfg.XFields{n});
                for iteration=1:datapoints
                    varyvars.(cfg.XFields{n})(:,:,m)= indepvars.(cfg.XFields{n})(:,:,m)*percentualrange(iteration);
                    eval(['[error output] = ' cfg.Model '(cfg,exvars,depvars,varyvars,r,1:130);']);
                    valuevector(iteration,z) = sum(output.TravelDemand(r,end,:))/1000;
                end
                labels{z} = [cfg.XFields{n} spacer usemodes{m}];
                z=z+1;
            end
        end
        
        globals.fid        = figure(4);set(4,'Units','centimeter')
        plot(percentualrange,valuevector);
        legend(labels);
        title(['Sensitivity analysis ' cfg.allregions{r}]);
        xlabel('Fractional change in variable');
        ylabel('Travel demand (thousand pkm/cap/yr)');
                pubgraph(globals.fid,10,1.5,'white');

        cfg.savename        = ['SA ' datestr(now,'yyyy-mm-dd HH-MM-SS' ) ' R ' cfg.allregions{r} ' Model ' cfg.Model];

    end
end

function savestuff(cfg)
    disp(['Saving graph as ' cfg.outputfolder cfg.OutputFile cfg.savename '.jpg']);
    saveas(gcf,[cfg.outputfolder cfg.OutputFile cfg.savename '.jpg']);
end